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Abstract. The notions of subgraph ccntrality and communicability, based on the exponential of 
the adjacency matrix of the underlying graph, have been effectively used in the analysis of undirected 
networks. In this paper we propose an extension of these measures to directed networks, and we apply 
them to the problem of ranking hubs and authorities. The extension is achieved by bipartization, i.e., 
the directed network is mapped onto a bipartite undirected network with twice as many nodes in order 
to obtain a network with a symmetric adjacency matrix. We explicitly determine the exponential of 
this adjacency matrix in terms of the adjacency matrix of the original, directed network, and we give 
an interpretation of centrality and communicability in this new context, leading to a technique for 
ranking hubs and authorities. The matrix exponential method for computing hubs and authorities 
is compared to the well known HITS algorithm, both on small artificial examples and on more 
realistic real- world networks. A few other ranking algorithms are also discussed and compared with 
our technique. The use of Gaussian quadrature rules for calculating hub and authority scores is 
discussed. 
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1. Introduction. In recent years, the study of networks has become central 
to many disciplines [5j HI HI [15j [16j ESI E3 EH] . Networks can be used to describe 
and analyze many different types of interactions, from those between people (social 
networks), to the flow of goods across an area (transportation networks), to links 
between websites (the WWW graph), and so forth. In general, a network is a set of 
objects (nodes) and the connections between them (edges). Often, research is focused 
on determining and describing important structural characteristics of a network or 
the interactions among its components. 

One common question in network analysis is to determine the most "important" 
nodes (or edges) in the network, also called node or vertex (edge) centrality. The 
interpretation of what is meant by "important" can change from application to appli- 
cation. Due to this, many different measures of centrality have been developed. For 
an overview, see [8] . A closely related notion is that of rank of a node in a network. 
There exist a number of definitions and algorithms for computing rankings; see, e.g., 
[H [3JJ [301 El El EJJ for up-to-date overviews. 

The main notion of node centrality considered in this paper, subgraph centrality, 
was introduced by Estrada and Rodriguez- Velazquez in . We refer readers to [5D] 
for the motivation behind this notion and for its name; see also the review article 
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[19] , and the discussion in section 2) The interpretation of ccntrality described in 
[19| applies mostly to undirected networks. However, many important real-world 
networks (the World Wide Web, the Internet, citation networks, food webs, certain 
social networks, etc.) are directed. One goal of this paper is to extend the notions 
of centrality and communicability described in [TTl [19] to directed networks, with an 
eye towards developing new ranking algorithms for, e.g., document collections, web 
pages, and so forth. We further compare our approach with some standard algorithms, 
such as HITS (see [55]) and a few others. Methods of quickly determining hub and 
authority rankings using Gauss-type quadrature rules are also discussed. 

2. Basic notions. Here we briefly review some basic graph-theoretic notions; 
we refer to |13j for a comprehensive treatment. A graph G = (V, E) is formed 
by a set of nodes (vertices) V and edges E formed by unordered pairs of vertices. 
Every network is naturally associated with a graph G = (V,E) where \V\ is the 
number of nodes in the network and E is the collection of edges between objects, 
E = | there is an edge between node i and node j}. The degree di of a vertex i 

is the number of edges incident to i. 

A directed graph, or digraph G = (V, E) is formed by a set of vertices V and 
edges E formed by ordered pairs of vertices. That is, € E ^> £ E. In the 
case of digraphs, which model directed networks, there are two types of degree. The 
in-degree of node % is given by the number of edges which point to i. The out-degree 
is given by the number of edges pointing out from i. 

A walk is a sequence of vertices v\, V2, ■ ■ ■ , Vk such that for 1 < i < k, there is an 
edge between Vi and Vi+± (or a directed edge from Vi to Vi+i in the case of a digraph). 
Vertices and edges may be repeated. A walk is closed if v± = Vk- A path is a walk 
consisting only of distinct vertices. 

A graph G is connected if every pair of vertices is linked by a path in G. A digraph 
is strongly connected if for any pair of vertices Vi and Vk there is a walk starting at Vi 
and ending at vj» . A digraph is weakly connected if the graph obtained by disregarding 
the orientation of its edges is connected. Unless otherwise specified, every digraph 
in this paper is simple (unweighted with no multiple edges or loops and connected). 
Note, however, that most of the techniques and results in the paper can be extended 
without difficulty to more general digraphs, in particular weighted ones. 

The adjacency matrix of a graph is a matrix A S ]Rl 1/ l >< l v 'l defined in the following 



Under the conditions imposed on G, A has zeros on the diagonal. If G is an 
undirected graph, A will be a symmetric matrix and the eigenvalues will be real. In 
the case of digraphs, A is not symmetric and may have complex (non-real) eigenvalues. 

3. Kleinberg's HITS algorithm. Here we briefly recall the classical Hypertext 
Induced Topics Search (HITS) algorithm, first introduced by J. Kleinberg in [29] . This 
algorithm provides the motivation for the extension of subgraph centrality to directed 
graphs given in section [5] 

3.1. The basic iteration. The HITS algorithm is based on the idea that in the 
World Wide Web, and indeed in all document collections which can be represented by 
directed networks, there are two types of important nodes: hubs and authorities. Hubs 
are nodes which point to many nodes of the type considered important. Authorities 
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1, if is an edge in G, 
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are these important nodes. From this comes a circular definition: good hubs are those 
which point to many good authorities and good authorities arc those pointed to by 
many good hubs. 

Thus, the HITS ranking relies on an iterative method converging to a stationary 
solution. Each node i in the network is assigned two non-negative weights, an authority 
weight x% and a hub weight yi. To begin with, each Xi and yi is given an arbitrary 
nonzero value. Then, the weights are updated in the following ways: 

xf = yf~ l) and x f ] for * = i»2,3... (3.1) 

The weights are then normalized so that Ylj( x j ) 2 = 1 an d YljiVj^) 2 = 

The above iterations occur sequentially and it can be shown that, under mild 
conditions, both sequences of vectors {x^} and {y^} converge as k — > oo. In 
practice, the iterative process is continued until there is no significant change between 
consecutive iterates. 

This iteration sequence shows the natural dependence relationship between hubs 
and authorities: if a node i points to many nodes with large x-valucs, it receives a 
large y-value and, if it is pointed to by many nodes with large y-valucs, it receives a 
large x- value. 

In terms of matrices, the equation (3.1) becomes: x^ k > — A T ?/ fe ~ 1 ) and — 
Ax^ , followed by normalization in the 2-norm. This iterative process can be ex- 
pressed as 

x^=c k A T Ax^ and yW = c' k AA T y^ , (3.2) 

where c k and c' k are normalization factors. A typical choice for the inizialization 
vectors x^°\ y(°> would be the constant vector 

see |21j . Hence, HITS is just an iterative power method to compute the dominant 
eigenvector for A T A and for AA T . The authority scores are determined by the entries 
of the dominat eigenvector of the matrix A T A, which is called the authority matrix 
and the hub scores are determined by the entries of the dominant eigenvector of 
AA T , called the hub matrix. Recall that the eigenvalues of both A T A and AA T are 
the squares of the singular values of A. Also, the eigenvectors of A T A are the right 
singular vectors of A, and the eigenvectors of AA T are the left singular vectors of A. 

3.2. HITS reformulation. In a digraph the adjacency matrix A is generally 
nonsymmetric, however, the two matrices used in the HITS algorithm (A T A and 
AA T ) are symmetric. Note that, setting 

a symmetric matrix is obtained. Now, 

, 2 _ / AA T \ , 3 _ / AA T A \ 
^-^0 A T A J ' V ATAAT / ' 

In general, 

A(A T A) k \ 

{A T A) k A T ) ■ 



A 2k = 



(AA T ) k 






(A T A) k 
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2k+l 



4 



M. Benzi, E. Estrada, and C. Klymko 



Applying HITS to this matrix A, A T = A so A T A = AA T = A 2 and introducing 
the vector = ( ^ J for k = 1, 2, 3, . . equation (3.2) becomes 

u^=AW-V = ( A f a t a )u^\ (3.3) 

followed by normalization of the two vector components of so that each has 2- 
norm equal to 1. Now, if A is an n X n matrix, A is 2n x 2n and vector is in 
R 2n . The first n entries of vS k ^ correspond to the hub rankings of the nodes, while 
the last n entries give the authority rankings. Under suitable assumptions (see the 
discussion in (32] Chapter 11.3]), as k — > oo the sequence {u^} converges to the 
dominant nonncgative eigenvector of A, which yields the desired hub and authority 
rankings. 

Hence, in HITS only information obtained from the dominant eigenvector of A 
is used. It is natural to expect that taking into account spectral information corre- 
sponding to the remaining eigenvalues and eigenvectors of A may lead to improved 
results. 

Among the limitations of HITS, we mention the possible dependence of the rank- 
ings on the choice of the initial vectors x^°\ y«, see [H] for examples of this, and 
the fact that HITS hub/authority rankings tend to be "degree-biased", i.e., they are 
strongly correlated with the out-/in-degrees of the corresponding nodes [14]. The 
latter property is in fact shared by most eigenvector-based rankings; for a discussion 
of this phenomenon in the case of scale-free graphs, see [35] . 

4. Subgraph centralities and communicabilities. In [19], the authors re- 
view several measures to rank the nodes in an undirected network A based on the use 
of matrix functions, such as the matrix exponential e A . The subgraph centrality |20j of 
node i is given by [e A ]u and the communicability |T7] between nodes i and j (i ^ j) is 
given by [e A ]ij. Nodes i corresponding to higher values of [e A ]u are considered more 
important than nodes corresponding to lower values. Large values of [e A ]ij indicate 
that information flows more easily between nodes i and j than between pairs of nodes 
corresponding to lower values of the same quantity. The Estrada index of the graph is 
given by Tr (e A ) = X)"=i [ eA \a- This index, which provides a global characterization 
of a network, is analogous to the partition function in statistical mechanics and plays 
an important role in the study of networks at the macroscopic level: quantities such 
as the natural connectivity, the total energy, the Helmholtz free energy and the entropy 
of a network can all be expressed in terms of the Estrada index |18] . 

Consider the power series expansion of e A , 

A 2 A 3 A k 

From graph theory, it is well known that if A is the adjacency matrix of an undirected 
graph, [A k ]ij = [A k ]ji counts the number of walks of length k between nodes i and 
j. Thus, the subgraph centrality of node i, [e A ]u, counts the total number of closed 
walks starting at node i, penalizing longer walks by scaling walks of length k by the 
factor jTj-. The communicability between nodes i and j, [e A ]ij, counts the number of 
walks between nodes i and j, again scaling walks of length A: by a factor of ^. 

It is worth mentioning that normalization of the diagonal entries of e A by Tr (e A ) 
yields a probability distribution on the nodes of the network, which can be given 
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the following interpretation: the ith diagonal entry of e^/Tr (e A ) is the probability 
of selecting any weighted self-returning (closed) walk that starts and ends at node i 
among all the weighted self-returning walks that start at any node and return to the 
same node. The weights used (factorial penalization) ensure that the shortest walks 
receive more weight than the longer ones: hence, the subgraph ccntrality of node i is 
proportional to the probability of finding a random walker walking "nearby" node i. 

Although the matrix exponential is certainly well-defined for any matrix, whether 
symmetric or not, the interpretation of the notion of subgraph ccntrality for directed 
networks can be problematic. To see this, consider the directed path graph consisting 
of n nodes, with edge set E = {(1, 2) , (2, 3) , . . . , (n — 1, n)} and adjacency matrix 
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The entries of e A are given by 



i/(i 
o, 



0! 



if 3 > i, 
else. 



In particular, the diagonal entries of e A are all equal to 1. Therefore, it is impos- 
sible to distinguish any of the nodes from the others on the basis of this centrality 
measure; yet, it is clear that the first and last node are rather special, and certainly 
more "peripheral" (less "central") than the other nodes. Also, we note that the prob- 
abilistic interpretation given above for undirected graphs is no longer meaningful for 
this example. Part of the problem, of course, is that the path digraph contains no 
closed walks. In the next section we show one way to extend the notion of subgraph 
centrality to digraphs that is immune from such shortcomings, and correctly differ- 
entiates between nodes in the example above. (On the other hand, it is interesting 
to note that the interpretation of the off-diagonal entries of e A in terms of commu- 
nicabilities is straightforward for the directed path. All entries of e A below the main 
diagonal are zero, reflecting the fact that information can only flow from a node to 
higher-numbered nodes. Also, the entries of e A decay rapidly away from the main 
diagonal, reflecting the fact that the "ease" of communication between a node and a 
higher numbered one decreases rapidly with the distance.) 

Another issue when extending the notions of subgraph centrality and communi- 
cability to directed graphs is that computational difficulties may arise. While the 
computations involved do not pose a problem for small networks, many real-world 
networks arc large enough that directly computing the exponential of the adjacency 
matrix is prohibitive. In [2], techniques for bounding and estimating individual entries 
of the matrix exponential using Gaussian quadrature rules are discussed; see also [6j 
and section |9] below. The ability to find upper and lower bounds for the entries re- 
quires that the matrix be symmetric, thus these bounds cannot be directly computed 
using the adjacency matrix of a directed network. Again, these difficulties can be 
circumvented using the approach discussed in the next section. 
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5. An extension to digraphs. Although the techniques described in [5] cannot 
be directly applied to non-symmetric matrices, setting 

> - ( jr t) ("I 

produces a symmetric matrix A and, thus, upper and lower bounds of individual en- 
tries of e" 4 can be computed. In Proposition [1] below we relate e -4 to the underlying 
hub and authority structure of the original digraph. By B> we denote the Moore- 
Penrose generalized inverse of matrix B. 



Proposition 1. Let A be as described in equation \5. 1}) . Then, 

1 cosh ( VXP") A (VWAj 1 sinh (vafa) 

^ sinh [y/WA^j (yWA) ^ A T cosh (VaTa) 



Proof. Let A = UT,V T be the SVD of the original (non-symmetric) adjacency 

A nm A ^ J i , ( U \ / S\ / U T 

matrix/!. 1 nen, A can be decomposed as A = \ q V / I E / \ V T 



Hence, 



Now, 



A I U \ ( S\ / U T 

6 = ( V CXP E ^ I' (5 - 2) 



S )= cosh ( S )+ sinh ( s 



cosh(E) 
cosh(E) 



Thus, 



exp 



E 
E 



sinh(E) 
sinh(E) 



cosh(E) sinh(E) 
sinh(E) cosh(E) 



(5.3) 



Putting together equations (|5.2|) and ()5.3 
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V 



cosh(E) sinh(E) 
sinh(E) cosh(E) 



U T 








cosh(%/Z4^) A(Vl^)^inh(\/I 5 "A) 



sinh Vl 5 ^ V^4 A T 



cosh ( VA^A 



\ 

The identities involving the off-diagonal blocks can be easily checked using the SVD 
of A. 
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5.1. Interpretation of diagonal entries. In the context of undirected net- 
works, the interpretation of the entries of the matrix exponential in terms of subgraph 
centralities and communicabilities is well-established, see e.g. [T^. In the case of di- 
rected networks and e A , things are not as clear. The network behind A can be thought 
of as follows: take the vertices from the original network A and make two copies of 
them, V and V. Then, undirected edges exist between the two sets based on the 
following rule: E' = \ there is a directed edge from i to j in the original network} . 

This creates a bipartite graph with 2n nodes: 1, 2, . . . , n, n + 1, n + 2, . . . , 2n. We 
denote by V(A) this set of nodes. The use of bipartization to treat rectangular and 
structurally unsymmetric matrices is of course standard in numerical linear algebra. 

In the undirected case, each node had only one role to play in the network: any 
information that came into the node could leave by any edge. In the directed case, 
there are two roles for each node: that of a hub and that of an authority. It is unlikely 
that a high ranking hub will also be a high ranking authority, but each node can still 
be seen as acting in both of these roles. In the network A, the two aspects of each 
node are separated. Nodes 1,2, ... ,n in V(A) represent the original nodes in their 
role as hubs and nodes n + 1, n + 2, . . . , 2n in l^(-A) represent the original nodes in 
their role as authorities. 

Given a directed network, an alternating walk of length k, starting with an out- 
edge, from node v\ to node Vk+i is a list of nodes v\, i>2, Vk+i such that there exists 
edge {vi,Vi+i) if i is odd and edge (vi+i,Vi) if i is even: 



Vi — > i>2 V3 — ► • " " 



An alternating walk of length k, starting with an in-edge, from node v\ to node Vk+i 
in a directed network is a list of nodes V\, V2, Vk+i such that there exists edge 
(i>i+i, Vi) if i is odd and edge (vi, if i is even: 



V\ <— V-2 — > t>3 • • • 



From graph theory (see also [TT| ) , it is known that [AA T A. . .]y (where there are k 
matrices being multiplied) counts the number of alternating walks of length k, starting 
with an out-edge, from node i to node j, whereas [A T AA T . . .]y (where there are k 
matrices being multiplied) counts the number of alternating walks of length k, starting 
with an in-edge, from node i to node j. That is, [(AA T ) k ]ij and \(A T A) k ]ij count the 
number of alternating walks of length 2k. 

In the original network A, if node i is a good hub, it will point to many good 
authorities, which will in turn be pointed at by many hubs. These hubs will also point 
to many authorities, which will again be pointed at by many other hubs. Thus, if i 
is a good hub, it will show up many times in the sets of hubs described above. That 
is, there should be many even length alternating walks, starting with an out-edge, 
from node i to itself. Giving a walk of length 2k a weight of -^yp these walks can be 
counted using the entry of the matrix 

AA T AA T AA T (AA T ) k 

/+ ^r + ^^ + -" + l2Aor + --- 

Letting A = UT,V T be the SVD of A, this becomes: 



S 2 S 4 S 2fc 

r+ 2! + 4T + - + (2*0T + - 



T 
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= Ucosh(E)U T = cosh(VAA T ) . 
The hub centrality of node i (in the original network) is thus given by 

[e A ]u = [cosh(VAAT)} H . 

This measures how well node i transmits information to the authoritative nodes in 
the network. 

Similarly, if node i is a good authority, there will be many even length alternating 
walks, starting with an in-edge, from node i to itself. Giving a walk of length 2k a 
weight of (2^)7 1 these walks can be counted using the (i, i) entry of cosh(V A T A). 

Hence, the authority centrality of node i is given by 

[e A ] n+ i tn+i = [cosh(V^ T A)] zl . 

It measures how well node i receives information from the hubs in the network. 

Note that the traces of the two diagonal blocks in e A are identical, so each accounts 
for half of the Estrada index of the bipartite graph. Also, recalling the well-known 
fact that the eigenvalues of A are icr^ where o~i denotes the singular values of A, we 
have 

n n n 

Tr (e A ) = e °* + Yl e ~ a ' = 2 cosh 

i—1 i—1 i—1 

an identity that can also be obtained directly from the expression for e A given in 
Proposition [TJ 

Returning to the example of the directed path graph with adjacency matrix A 
given by (|4.2p . one finds that using the diagonal entries of e A to rank the nodes 
gives node 1 as the least authoritative node, and node n as the one with the lowest 
hub ranking, with all the other nodes being tied. Thus we see that, while e A fails 
to differentiate between the nodes of this graph, using e A yields a very reasonable 
hub/authority ranking of the nodes. 

5.2. Interpretation of off-diagonal entries. Although not used in the re- 
mainder of this paper, for the sake of completeness we give here an interpretation of 
the off-diagonal entries of e A . As we will see, this interpretation is rather different 
from the one usually given for the off-diagonal entries of e A , and provides information 
of a different nature on the structure of the underlying graph. 

In discussing the off-diagonal entries of A, there are three blocks to consider. 
First, there are the off-diagonal entries of the upper- left block, cos h(v / H ? ), then 
there are the off-diagonal entries of the lower-right block, cos h(V A T A). Finally, there 

is the off-diagonal block, A A T Aj sinh (V A T Aj (the fourth block in e A being its 
transpose). 

From section [O] [e A ]ij = [cosh(V A A T )]ij , 1 < i, j < n, counts the number 
of even length alternating walks, starting with an out-edge, from node i to node j, 
weighting walks of length 2k by a factor of j^y - When i ^ j, these entries measure 
how similar nodes i and j are as hubs. That is, if nodes i and j point to many of the 
same nodes, there will be many short even length alternating walks between them. 

The hub communicability between nodes i and j, 1 < i,j < n, is given by 



[e-% = [cosh(v / Z4 ? )]. i 
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This measures how similar nodes i and j are in their roles as hubs. That is, a larger 
value of hub communicability between nodes i and j indicates that they point to many 
of the same authorities. In other words, they point to nodes which are authorities on 
the same subjects. 

Similarly, [e A ] n +i, n +j = [coah(V A T A)]ij , 1 < i,j < n, counts the number of 
even length alternating walks, starting with an in-edge, from node i to node j, also 
weighing walks of length 2k by a factor of jjj^r- When i ^ j, these entries measure 
how similar the two nodes are as authorities. If i and j are pointed at by many of the 
same hubs, there will be many short even length alternating walks between them. 

The authority communicability between nodes i and j, 1 < i,j, < n, is given by 

[e A ] i+nd+n = [cosh{V A T A)] tj 

This measures how similar nodes i and j are in their roles as authorities. That is, a 
larger value of authority communicability between nodes i and j means that they are 
pointed to by many of the same hubs and, as such, are likely to contain information 
on the same subjects. 

Let us now consider the off-diagonal blocks of A. Here, [sinh(-\/ A T A)]ij counts 
the number of odd length alternating walks, starting with an out-edge, from node i to 
node j, weighing walks of length 2k +1 by ( 2 k+i)\ • This measures the communicability 
between node i as a hub and node j as an authority. 

The hub- authority communicability between nodes i and j (that is, the communi- 
cability between node i as a hub and node j as an authority) is given by: 

[e A ] i>n+j = [A (VaTa)* sinh (v^l)]^ 



= [sinh (VI 5 ^) (71^ A T ]ji = [e A ] n+h i. 

A large hub-authority communicability between nodes i and j means that they are 
likely in the same "part" of the directed network: node i tends to point to nodes that 
contain information similar to that on which node j is an authority. 

5.3. Relationship with HITS. As described in 13.21 the HITS ranking of nodes 
as hubs and authorities uses only information from the dominant eigenvector of A. 
Here we show that when using the diagonal of e A , we exploit information contained 
in all the eigenvectors of A] moreover, the HITS rankings can be regarded as an 
approximation of those given by the diagonal entries of e A . 

Assume the eigenvalues of A can be ordered as Ai > A 2 > A3 > • • • > \2 n - 
Then, A can be written as A = X)i=i ^i u i u f where m, 1*2, • • • , ui n are the normalized 
eigenvectors of A. Taking the exponential of A, we get: 

2n 2n 

e A = e Xi Ujuf = e Xl uiu{ + e x 'Uju[. 
Now, the hub and authority rankings come from the diagonal entries of e A : 

2n 

diag (e A ) = e Al diag (uiuj) + e A, diag (wjwf ). 

i=2 
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Rcscaling the hub and authority scores by e 1 does not alter the rankings; hence, we 
can instead consider 

2n 

diag (e~ Xl e A ) = diag (e" A " Al/ ) = diag (umf) + ^ e A *- Al diag (u t uf). 

i=2 

Now, the diagonal entries of the rank-one matrix u\u\ are just the squares of the 
(nonnegativc) entries of the dominant eigenvector of A; hence, the rankings provided 
by the first term in the expansion of e A in the eigenbasis of A are precisely those 
given by HITS. 

It is also clear that if Ai ^> A2, then the rankings provided by the diagonal entries 
of e A are unlikely to differ much from those of HITS, since the weights e Ai ~ Al will 
be tiny, for all i = 2, . . . , 2n. Conversely, if the gap between Ai and the rest of the 
spectrum is small (Xi ~ A2), then the contribution from the remaining eigenvectors, 
yij—v e Ai_Al diag (utuf), may be non- negligible relative to the first term and therefore 
the resulting rankings could differ significantly from those obtained using HITS. In 
section [5] we will see examples of real networks illustrating both scenarios. 

Summarizing, use of the matrix exponential for ranking hubs and authorities 
amounts to using the (squared) entries of all the eigenvectors of .4, weighted by the 
exponential of the corresponding eigenvalues. Of course, in place of the exponential, 
a number of other functions could be used; see the discussion in the next section. 
Although using an exponential weighting scheme may at first sight appear to be 
arbitrary, its use can be rigorously justified; see the discussion in the next section, and 
[T8] for a thorough treatment in the context of undirected graphs. As shown above, the 
HITS ranking scheme uses the leading term only, corresponding to the approximation 
p A e ^i UlU J _ Between these two extremes one could also use approximations of the 
form 

fc 

e A *J2 eKuiU i> ( 5 - 4 ) 
i=i 

where 1 < k < n; indeed, in most cases of practical interest a modest value of k (<g; n) 
will be sufficient for a very good approximation, since the eigenvalues of A are often 
observed to decay rapidly from a certain index k onward. We return on this topic in 
section [5] 

6. Other ranking schemes. In this section we discuss a few other schemes that 
have been proposed in the literature, and compare them with the hub and authority 
centrality measures based on the exponential of A. 

6.1. Resolvent-based measures. Besides the matrix exponential, another func- 
tion that has been successfully used to define centrality and communicability measures 
for an undirected network is the matrix resolvent, which can be defined as 

R(A; c) = (/ - cA)- 1 =I + cA + c 2 A 2 + ■■■+ c k A k + • • ■ , 

with < c < l/A max (A). This approach was pioneered early on by Katz [28], and 
variants thereof have since been used by numerous authors; see, e.g., [6j [18j HU 
[23l |4T] . Here A is the symmetric adjacency matrix of the undirected network. The 
condition on the parameter c ensures that R(A; c) is well defined (i.e., that / — cA is 
invertible and the geometric series converges to its inverse) and nonnegative; indeed, 
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7 — cA will be a nonsingular ill-matrix. It is hardly necessary to mention the close 
relationship existing between the resolvent and the exponential function, which can 
be expressed via the Laplace transform. For the adjacency matrix A of a bipartite 
graph given by (|5.ip . the resolvent is easily determined to be 



The condition on c can be expressed as < c < l/oi, where o\ = \A\i denotes 
the largest singular value of A, the adjacency matrix of the undirected network. This 
ensures that the matrix in (|6.1[) is well-defined and nonncgative, with positive diagonal 
entries. The diagonal entries of (7 — c 2 'AA T )~ 1 provide the hub scores, those of 
(7 — c 2 A T A)~ 1 the authority scores. A drawback of this approach is the need to 
select the parameter c, and the fact that different values of c may lead to different 
rankings. We have performed numerical experiments with this approach and we 
found that for certain values of c, particularly those close to the upper limit l/oi, 
the hub and authority rankings obtained with the resolvent function are not too 
different from those obtained with the matrix exponential. However, not surprisingly, 
as the value of c is reduced, one obtains hub and authority rankings that are strongly 
correlated with the out- and in-degree of the nodes, respectively^!] Overall, because 
the resolvent tends to weigh short walks more heavily than the exponential, and 
since longer walks contribute relatively little to the centrality scores, it is fair to 
say that the exponential is less "degree biased" than the resolvent function. Also, 
since the exponential rankings do not depend on a tuneable parameter, they provide 
unambiguous rankings. 

We note that "Katz" authority and hub scores may also be obtained by consid- 
ering the column and row sums of the (nonsymmetric) matrix resolvent (7 — cA)^ 1 , 
where A is the adjacency matrix of the original digraph and c > is again assumed 
to be small enough for the corresponding Neumann series to converge. Indeed, the 
row sums of (7 — cA)" 1 count the number of (weighted) walks out of each node, while 
the column sums count the number of (weighted) walks into each node. Denoting by 
1 the vector of all ones, hub and authority rankings can be obtained by solving the 
two linear systems 



respectively. Here the parameter c must satisfy < c < l/p(A), where p(A) denotes 
the spectral radius of A. The results of numerical experiments comparing the Katz 
scores with those based on the exponential of A are given in section [8J Here we 
observe that these Katz scores are also dependent on the choice of the parameter c, 
and similar considerations to those made for (7 — cA) -1 apply. 

A natural analogue to this approach is the use of row and column sums of the ex- 
ponential e A to rank hubs and authorities. Some results obtained with this approach 
are discussed in section [8j We note that this method is different from the Exponenti- 
ated Inputs HITS Method of [21] . The latter method is a modification to HITS which 
was developed in order to correct the issue of non unique results in certain networks. 
If the dominant eigenvalue of A T A (and, consequently, of AA T ) is not simple, then 
the corresponding eigenspace is multidimensional. This means that the choice of the 




(6.1) 



(7 - cA)y = 1 and (7 - cA 



(6.2) 



Note that if c is taken too small, then the resolvent approaches the identity matrix and it 
becomes impossible to have meaningful rankings of the nodes. 
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initial vector affects the convergence of the HITS algorithm and different hub and 
authority vectors can be produced using different initial vectors. This can occur only 
when A T A is reducible, that is, when the original network is not strongly connected. 
In |21j . Farahat et al. propose a modification to the HITS algorithm which amounts 
to replacing A and A T with e A — I and (e A — I) T in the HITS iteration. They note 
that, as long as the original network is weakly connected, the dominant eigenvalue 
of (e A — I) T (e A — I) is simple. Thus, HITS with this exponentiated input produces 
unique hub and authority rankings. However, a result of this replacement is that 
nodes with zero in-degree (or a low in-degree) are less important in the calculation of 
authority scores than nodes with a high in-degree. When there are many nodes with 
zero in-degree or whose edges point to only a few other nodes, dropping these edges 
can greatly affect the HITS rankings. An obvious disadvantage of this algorithm is 
its cost, since it requires iterating with a matrix exponential and its transpose. It 
can be implemented using only matrix- vector products involving A and A T by means 
of techniques, like Krylov subspacc methods, for evaluating the action of a matrix 
function on a given vector; see, e.g., [26l Chapter 13]. This approach leads to a nested 
iteration scheme, with HITS as the outer iteration and the Krylov method as the 
inner one. Generally speaking, we have found HITS with exponentiated inputs to be 
less reliable and more expensive than the other methods considered in this paper. We 
refer to [3] for additional discussion and some examples. 

6.2. PageRank and Reverse PageRank. As is well known, the (now) clas- 
sical PageRank algorithm provides a means of finding the authoritative nodes in a 
digraph. In PageRank, the importance of a node v is determined by the importance 
of the nodes pointing to v. In the most basic formulation, the rank of v is given by 

(«) 



it. 

ueB v 1 1 

where B v = {u : there is a directed link from u to v} and \u\ is the out-degree of 
u. The ranks of the nodes are computed by initially setting, say, r(v) = ^ (where n 
is the size of the network) and iteratively computing the rankings until convergence. 
This can also be written as 

TT^TT^P, k = l,2,... (6.4) 

where -Kk is the vector of node ranks at iteration k and P is the matrix given by 

J 1/| Vi | , if there is a directed edge from Vi to Vj , 
Pij ~ { 0, else. 

Here, P can be viewed as a probability transition matrix, where is the probability 
of traveling from node Vi to node Vj along an edge and the iterations can be understood 
as the evolution of a Markov chain modeling a random walk on the network. 

However, for an arbitrary network, there is no guarantee that the PageRank algo- 
rithm will converge. If there are nodes with zero out-degree, P will not be stochastic. 
To correct this, the matrix P is used, where each zero row of P is replaced with e T /n. 
Although this guarantees that the algorithm will converge, it does not guarantee the 
existence of a unique solution. Even with the augmentation, P might still be a re- 
ducible matrix, corresponding to a reducible Markov chain. When this happens, there 
are rank sinks, i.e., nodes in which the random walk will become trapped and, subse- 
quently, these nodes will receive a disproportionately high rank. However, if P were 
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irreducible, there would be no rank sinks and the Perron-Frobenius theorem would 
guarantee that the Markov chain had a unique, positive stationary distribution. 

The standard way to form a stochastic, irreducible PagcRank matrix P is to 
introduce the rank-1 matrix E = ee T /n and to consider instead of P the convex 
combination 



where a is a constant with < a < 1. The coefficient 1 — a is a measure of the 
tendency of a person surfing the web to jump from one page to another without 
following links. In practice, a frequently recommended value is a = 0.85. For a more 
comprehensive overview of the PageRank algorithm, see [23l [27l |3TJ [32] . 

It was pointed out in [22] that applying PageRank to the digraph obtained by 
reversing the direction of the edges provides a natural way to rank the hubs; this 
is usually referred to as Reverse PageRank. In other words, authority rankings are 
obtained by applying PagcRank to the "Google" matrix derived from A, and hub 
rankings are obtained by the same process applied to A T . Like HITS, PageRank and 
Reverse PageRank are eigenvector-based ranking algorithms that do not take into ac- 
count information about the network contained in the non-dominant eigenvectors. As 
already mentioned, it has been argued |35j that eigenvector-based algorithms tend to 
be degree-biased. Furthermore, like the Katz-type algorithms, the rankings obtained 
depend on the choice of a tuneable damping parameter. While the success of PageR- 
ank in finding authoritative nodes is well known and very well documented, much 
less is known about the effectiveness of Reverse PageRank in identifying hubs; some 
references are [TJ [101 E2J E3] . We present the results of a few numerical experiments 
with PageRank and Reverse PagcRank in section [5] 

7. Examples. In this section and the next we illustrate the proposed method on 
some simple networks of small size, as well as on some larger data sets corresponding to 
real networks. We also compare our approach with HITS and other rankings schemes, 
including Katz, PageRank and Reverse PagcRank. 

7.1. Small digraphs. In this section we compare out and in-degree counts, 
HITS, and our proposed method to obtain hub and authority rankings in a few small 
digraphs. The purpose of this section is mostly pedagogical. 

7.1.1. Example 1. Consider the small directed network in Fig. 17.11 (left panel). 
The adjacency matrix is given by 



The corresponding bipartite graph is shown in Fig. 17.11 (right panel). If hubs and 
authorities are determined simply using in-degree and out-degree counts, the result is 



P = aP + (1 -a)E, 



(6-5) 



A = 



( 1 1 \ 
10 10 
10 1 

\ 1 J 



as follows: 



node 



out-degree 



in-degree 



1 

2 
3 
4 



2 
2 
2 
1 



1 

3 
2 
1 
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Fig. 7.1. The original directed network from Example 1, with adjacency matrix A (left) and 
the bipartite network with adjacency matrix A (right). 

Under this ranking, the hub ranking of the nodes is: {1, 2, 3 (tie); 4}. The au- 
thority ranking of the nodes is: {2; 3; 1, 4 (tie)}. We obtain somewhat different results 
using the HITS algorithm. The eigenvectors of AA T and A T A corresponding to the 
largest eigenvalue A max ~ 3.9563, which is simple, yield the following rankings for 
hubs and authorities: 



node 


hub rank 


authority rank 


1 


.3383 


.0965 


2 


.1729 


.4618 


3 


.2798 


.2854 


4 


.2091 


.1562 



Here, the ranking for hubs is: {1; 3; 4; 2}. The ranking for authorities is: {2; 3; 4; 1}. 
Note that node 2, which was given a top hub score by looking just at the out-degrees, 
is judged by HITS as the node with the lowest hub score. 

Using e" 4 as described above, the rankings for hub centralities and authority 
centralities are: 



node 


hub ccntrality = [e A ]u 


authority centrality = [c^+i^+j 


1 


2.3319 


1.5906 


2 


2.2289 


3.0209 


3 


2.2812 


2.2796 


4 


1.6414 


1.5922 



With this method, the hub ranking of the nodes is: {1;3;2;4}. The authority 
ranking is: {2;3;4;1}. On this example, our method produces the same authority 
ranking as HITS. The hub ranking, however, is slightly different: both methods iden- 
tify node 1 as the one with the highest hub score, followed by node 3; however, our 
method assigns the lowest hub score to node 4 rather than node 2. This is arguably 
a more meaningful ranking. 

7.1.2. Example 2. Consider the small directed network in Fig. 17.21 (left panel). 
The adjacency matrix is given by 

/ 1 \ 

10 1 

10 0' 
\ 1 J 
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Fig. 7.2. The original directed network from Example 2, with adjacency matrix A (left) and 
the bipartite network with adjacency matrix A ( right ) . 

The corresponding bipartite graph is shown in Fig. 17.21 (right panel). If hubs 
and authorities are determined only using in-degrees and out-degrees, the result is as 
follows: 



node 


out-degree 


in-dcgrcc 


1 


1 


1 


2 


2 


2 


3 


1 


1 


4 


1 


1 



Under this criterion, the hub and authority rankings are both {2; 1, 3, 4 (tie)}. 
While it is intuitive that node 2 should be given a high score (both as an authority 
and as a hub), just looking at the degrees does not allow one to distinguish the 
remaining nodes. 

Consider now the use of HITS. The largest eigenvalue of AA T (and A T A) is 
Amax = 2 and it has multiplicity two. Thus, different starting vectors for the HITS al- 
gorithm may produce different rankings, as discussed in [21j . Starting from a constant 
authority vector x^°\ as suggested in [55], produces the following scores: 



node 


hub rank 


authority rank 


1 


.0000 


.3333 


2 


.5000 


.3333 


3 


.2500 


.0000 


4 


.2500 


.3333 



The ranking for hubs is: {2; 3, 4 (tie); 1}. The ranking for authorities is the following: 
{1,2, 4 (tie); 3}. 

If the ranking is determined using e A as described above, the resulting scores are: 



node 


hub ccntrality = [e A ]u 


authority centrality = [e-^+j^+j 


1 


1.5431 


1.5891 


2 


2.1782 


2.1782 


3 


1.5891 


1.5431 


4 


1.5891 


1.5891 



With this method, the hub ranking of the nodes is the same as in HITS: {2; 3, 4 (tie); 1}. 
However, in the authority ranking, node 2 is the clear winner rather than being part 
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Fig. 7.3. The original directed network from Example 3, with adjacency matrix A (left) and 
the bipartite network with adjacency matrix A (right). 

of a three-way tie for first place: {2; 1,4 (tie); 3}. In this example, the method based 
on the matrix exponential is able to identify a top authority node by making use of 
additional spectral information. 

7.2. Example 3. Let G be the small directed network in Fig. 17. 31 The adjacency 
matrix is given by 

/000000\ 
1 
._ 10 

1 ■ 
1 
\ 1 1 1 1 / 

If hubs and authorities are determined using only in-degrees and out-degrees, the 
result is: 



node 


out-degree 


in-degree 


1 





4 


2 


1 


1 


3 


1 


1 


4 


1 


1 


5 


1 


1 


6 


4 






The hub ranking of the nodes using degrees is: {6; 2,3,4,5 (tie); 1}. The authority 
ranking is {1; 2,3,4,5 (tie); 6}. 

If the HITS algorithm is used, the resulting rankings are similar, but not exactly 
the same. Starting with a constant authority vector , the results are: 



node 


hub rank 


authority rank 


1 


.000 


.200 


2 


.125 


.200 


3 


.125 


.200 


4 


.125 


.200 


5 


.125 


.200 


6 


.500 


.000 



The hub ranking of the nodes is: {6; 2, 3, 4, 5 (tie); 1}. The authority ranking is: 
{1,2,3,4,5 (tie); 6}. Here, HITS does not differentiate between node 1 and nodes 2, 3, 
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4, and 5 in terms of the authority score, even though node 1 has by far the highest 
in-degree. This appears as a failure of HITS, since it is intuitive that node 1 should 
be regarded as very authoritative. 

When e A is used to calculate the hub and authority scores, node 1 does get a 
higher authority ranking than all the other nodes: 



node 


hub centrality = [e A ]u 


authority centrality = [e-^g+i^+i 


1 


1.0000 


3.7622 


2 


1.6905 


1.6905 


3 


1.6905 


1.6905 


4 


1.6905 


1.6905 


5 


1.6905 


1.6905 


6 


3.7622 


1.0000 



Note that, if desired, the value 1 can be subtracted from these scores since it does 
not affect the relative ranking of the nodes. The hub ranking is {6; 2,3,4,5 (tie); 1}, 
and the authority ranking is: {1; 2,3,4,5 (tic); 6}. 

8. Application to web graphs. Similarly to HITS, and in analogy to subgraph 
centrality for undirected networks, the rankings produced by the values on the diago- 
nal of e A can be used to rank websites as hubs and authorities in web searches (many 
other applications are of course also possible) . Three of the data sets considered here 
are small web graphs consisting of web sites on various topics and can be found at 
[4"0] along with the website associated with each node; see also [7]. The experiments 
for this paper were run on the "Expanded" version of the data sets. Each data set is 
named after the corresponding topicH In addition, wc include results for the wb-cs- 
stanford data set from the University of Florida Sparse Matrix Collection [12] . This 
digraph represents a subset of the Stanford University web. In this section, the hub 
and authority rankings obtained from e A are compared with those from HITS, Katz 
(using (|6.2|) with c = l/(p(A) + 0.1)), the row and column sums of the exponential 
e A of the nonsymmetric matrix A, and PageRank/Reverse PageRank. For the latter 
we use the standard value a = 0.85 for the damping parameter. All experiments are 
performed using Matlab Version 7.9.0 (R2009b) on a MacBook Pro running OS X 
Version 10.6.8, a 2.4 GHZ Intel Core i5 processor and 4 GB of RAM. For the purpose 
of these tests we use the built-in Matlab function expm to compute the matrix expo- 
nentials, and backslash to compute the Katz scores. Other approximations of e A are 
discussed in section [9] 

8.1. Abortion data set. The abortion data set contains n = 2293 nodes and 
m = 9644 directed edges. The expanded matrix A = ( ^ T ^ ~\ has order N = 

2n = 4586 and contains 2m = 19288 nonzeros. The maximum eigenvalue of A is 
\ N ps 31.91 and the second largest eigenvalue is Ajv-i ~ 26.04. In this matrix, the 
largest eigenvalue is fairly well-separated from the second largest so that one would 
expect the HITS rankings (which only use information from the dominant eigenpair 
of A) to be reasonably close to the rankings from e A (which use information from 
all of the eigenvalues and corresponding eigenvectors). A plot of the eigenvalues of 
the expanded abortion data set matrix can be found in Fig. 18.11 Note the high 



2 It should be noted, however, that in the node list for the adjacency matrix, the node labeling 
begins with 1 and in the list of websites associated with the nodes found at |40| , node labeling begins 
at 0. Thus, node i in the adjacency matrix is associated with website i — 1. 
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500 1000 1500 2000 2500 3000 3500 4000 4500 5000 
Fig. 8.1. Plot of the eigenvalues of the expanded abortion matrix A. 



Table 8.1 

Top 10 hubs of the abortion web graph, ranked using \e^\u, HITS, Katz, e A row sums and 
Reverse PageRank with a = 0.85. 



\p ] ii 


HITS 


Katz 


4 

e rs 


RPR 


48 


48 


80 


80 


125 


1021 


1006 


1431 


1431 


2184 


1007 


1007 


1432 


1432 


79 


1006 


1021 


1387 


1426 


81 


1053 


1053 


1388 


1425 


48 


1020 


1020 


1389 


1415 


1424 


987 


960 


1397 


1388 


1447 


990 


968 


1425 


1389 


78 


985 


969 


1426 


1397 


134 


989 


970 


1415 


1387 


1445 



multiplicity of the zero eigenvalue in this matrix, as well as in the adjacency matrices 
of the computational complexity and death penalty data sets. Also, quite a few of the 
nonzero eigenvalues are rather small. Due to this, the numerical rank of the matrix 
is very low, a property that can be exploited when estimating the entries of e A using 
Lanczos-based methods; see section [9] for further discussion on this. 

The top 10 hubs and authorities of the abortion data set, as determined using 
the diagonal entries of e* 4 , HITS with constant initial vector, the row/column sums of 
(7 — cA)^ 1 ("Katz"), the row/column sums of e A and Reverse PageRank/PageRank 
are shown in Tables 18.11 and 18.21 We observe that there is a good deal of agreement 
between the e A rankings and the HITS ones: indeed, both methods identify the 
websites labeled 48, 1021, 1007, 1006, 1053, 1020 as the top 6 hubs, and both pick 
web site 48 as the top one. Also, there are 7 web sites identified by both methods as 
being among the top 10 authorities. The top authority identified by HITS is ranked 
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Table 8.2 

Top 10 authorities of the abortion web graph, ranked using [e^\a, HITS, Katz, e A column sums 
and PageRank with a = 0.85. 



r Ai 

[e A )ii 


HITS 


Katz 


e cs 


PR 


967 


939 


1430 


1430 


1609 


958 


958 


1387 


1387 


1941 


939 


967 


1425 


1425 


1948 


962 


961 


1426 


1426 


1608 


963 


962 


1429 


1417 


587 


964 


963 


1396 


1409 


1610 


961 


964 


1405 


1429 


2045 


965 


965 


1406 


1406 


317 


966 


966 


1409 


1396 


2191 


587 


1582 


1417 


1405 


753 




200 400 600 800 1000 1200 1400 1600 1800 

Fig. 8.2. Plot of the eigenvalues of the expanded computational complexity matrix A. 



third by e , and conversely the top authority identified by e A is third in the HITS 
ranking. The Katz rankings and those based on e A show considerable agreement 
with one another, but are very different from the HITS ones and from those based 
on e A . Node 48, which is the top-ranked hub according to HITS and e A , is now 
not even among the top 100. Conversely, node 80, which is ranked the top hub by 
Katz and e A , is not in the top 100 nodes according to HITS or to e . This is not 
too surprising, since the metrics based on A and those based on A are obtained 
by counting rather different types of graph walks. Finally, for this network Reverse 
PageRank and PageRank return rankings with almost no overlap with any of the 
other methods. 

8.2. Computational complexity data set. The computational complexity 
data set contains n = 884 nodes and m = 1616 directed edges. The expanded matrix 
A has order N = 2n = 1768 and contains 2m = 2232 nonzcros. The maximum eigen- 
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Table 8.3 

Top 10 hubs of the computational complexity web graph, ranked using [e A ]a, HITS, Katz, e A 
row sums and Reverse PageRank with a = 0.85. 



r Ai 

[er]u 


HITS 


Katz 


e rs 


RPR 


57 


57 


56 


57 


57 


17 


634 


709 


56 


56 


644 


644 


57 


17 


17 


643 


721 


697 


51 


51 


634 


643 


705 


634 


21 


106 


544 


690 


21 


11 


119 


632 


714 


255 


255 


529 


801 


708 


173 


12 


86 


640 


712 


709 


13 


162 


639 


715 


45 


45 



value of A is An ~ 10.93 and the second largest eigenvalue is Ajv-i ~ 9.86. Here, the 
(relative) spectral gap between the first and the second eigenvalue is smaller than in 
the previous example; consequently, we expect the rankings produced using e A and 
HITS to be less similar than for the abortion data set. A plot of the eigenvalues of 
the expanded computational complexity data set matrix can be found in Fig. 18.21 

The top 10 hubs and authorities of the computational complexity data set, de- 
termined by the various ranking methods, can be found in Tables 18.31 and 18.41 As 
expected, we see less agreement between HITS and the diagonals of e A . Concerning 
the hubs, both methods agree that the web site labelled 57 is by far the most impor- 
tant hub on the topic of computational complexity. However, the method based on 
e A identifies as the second most important hub the web site corresponding to node 17, 
which is ranked only 39th by HITS. The two methods agree on the next three hubs, 
but after that they return completely different results. The difference is even more 
pronounced for the authority rankings. The method based on e A clearly identifies 
web site 1 as the most authoritative one, whereas HITS relegates this node to 8th 
place. The top authority acording to HITS, web site 719, places 5th in the ranking 
obtained by e A . The two methods agree on only two other web sites as being in the 
top 10 authorities (717 and 727). The Katz rankings and those based on e A show 
little overlap for this data set, although node 57 is clearly considered an important 
hub by all measures. A natural question is how much these results are affected by the 
choice of the parameter c used to compute the Katz scores. We found experimentally 
that, in contrast to the situation for the other data sets, small changes in the value of 
c can significantly affect the Katz ranking for this particular data set. Changing the 
value of c to c = l/(p(A) + 0.3) results in hub and authority rankings that are much 
closer to those given by the column/row sums of e A . The potential sensitivity to c is 
a clear drawback of the Katz-based approach compared to the methods based on the 
matrix exponential. Coming to (Reverse) PageRank, it is interesting to note that for 
this data set it provides rankings that are at least in partial agreement with some of 
the other measures, especially those based on e A . Looking at the authority scores, 
we also notice a good degree of overlap among all methods, except HITS. Due to the 
small spectral gap, HITS is probably the least reliable of all ranking methods on this 
particular data set. 
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Table 8.4 

Top 10 authorities of the computational complexity web graph, ranked using [e J 
e A column sums and PageRank with a = 0.85. 



HITS, Katz, 



r Ai 


HITS 


Katz 


e A cs 


PR 


1 


719 


688 


673 


673 


315 


717 


685 


1 


664 


673 


727 


673 


664 


534 


148 


723 


690 


534 


45 


719 


808 


56 


45 


2 


717 


735 


686 


473 


1 


2 


737 


664 


315 


376 


45 


1 


1 


376 


341 


727 


722 


45 


688 


50 


534 


770 


534 


599 


51 




-10 
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Fig. 8.3. Plot of the eigenvalues of the expanded death penalty matrix A. 



8.3. Death penalty data set. The death penalty data set contains n = 1850 
and m = 7363 directed edges. The expanded matrix A has order N = 2n = 3700 and 
contains m = 14726 nonzeros. The maximum eigenvalue of A is A^r ~ 28.02 and the 
second largest eigenvalue Aat_i ps 17.68. In this case, the largest and second largest 
eigenvalues are quite far apart, and the relative gap is larger than in the previous 
examples. A plot of the eigenvalues of the expanded death penalty matrix can be 
found in Fig. IO 

Due to the presence of a large spectral gap, much of the information used in 
forming the rankings of is also used in the HITS ranking, and we expect the two 
methods to produce similar results; see section 15.31 Indeed, as shown in Table 18.51 
(hubs) and Table (authorities), in this case the top 10 rankings produced by the 
two methods are actually identical. 

Looking at the Katz scores and those based on e A , we see in this case a great 
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Table 8.5 

Top 10 hubs of the death penalty web graph, ranked using [eA\u, HITS, Katz, e A row sums and 
Reverse PageRank with a = 0.85. 



r Ai 

[er]u 


HITS 


Katz 


e rs 


RPR 


210 


210 


1632 


1632 


210 


637 


637 


133 


133 


1632 


413 


413 


1671 


1671 


70 


1586 


1586 


552 


552 


95 


552 


552 


1651 


1651 


135 


462 


462 


1673 


210 


133 


930 


930 


1328 


1673 


55 


542 


542 


1653 


1653 


958 


618 


618 


210 


1328 


1077 


1275 


1275 


1709 


1709 


315 



Table 8.6 

Top 10 authorities of the death penalty web graph, ranked using \e^\a, HITS, Katz, e A column 
sums and PageRank with a = 0.85. 



[e A h 


HITS 


Katz 


e cs 


PR 


4 


4 


1632 


1632 


993 


1 


1 


1662 


1662 


667 


6 


6 


1697 


1697 


3 


7 


7 


1689 


1689 


736 


10 


10 


1653 


1653 


735 


16 


16 


1671 


1671 


1632 


2 


2 


1675 


1675 


42 


3 


3 


1684 


1684 


1 


44 


44 


798 


789 


4 


27 


27 


1652 


1654 


1212 



deal of overlap between these two, but almost completely different rankings compared 
to HITS and e A (although node 210 is clearly an important hub by any standard). 
Note that node 1632 is both the top hub and the top authority according to Katz 
and to e A . PageRank and Reverse PageRank show a limited amount of overlap with 
the other measures; nevertheless, nodes 210 and 1632 are also found to be important 
hubs and nodes 1632, 1 and 4 are found to be authoritative, in agreeemnt with some 
of the other measures. 

8.4. Stanford web graph. The wb-cs-stanford data set from the University of 
Florida sparse matrix collection contains n = 9914 nodes and m = 36854 directed 
edges. The expanded matrix A has order TV = 2n = 19828 and contains m = 73708 
nonzeros. The maximum eigenvalue of A is Ajv ~ 38.38 and the second largest is 
Aat_i rj 32.12, hence there is a sizeable gap. Tables [877118 . 81 report the results obtained 
with the various ranking schemes. 

The first thing to observe is the remarkable agreement between the HITS, e A , 
Katz, and e A rankings of both hubs and authorities. This in stark contrast with the 
results for the previous three data sets. Moreover, many of the nodes that are ranked 
highly as hubs are also ranked highly as authorities. A plausible explanation of these 
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Table 8.7 

Top 10 hubs of the wb-cs-stanford web graph, ranked using [e A ]u, HITS, Katz, e A row sums 
and Reverse PageRank with ct = 0.85. 



r Ai 

[er]u 


HITS 


Katz 


e rs 


RPR 


6562 


6562 


6562 


6562 


251 


6838 


6838 


6837 


6837 


252 


6840 


6837 


6838 


6838 


253 


6837 


6839 


6839 


6839 


254 


6839 


6840 


6840 


6840 


271 


6616 


6616 


6669 


6669 


2240 


6765 


6615 


6668 


6668 


2241 


6615 


6765 


6670 


6670 


2242 


6669 


6669 


6616 


6616 


2243 


6731 


6731 


6615 


6615 


348 



Table 8.8 

Top 10 authorities of the wb-cs-stanford web graph, ranked using [e^]«, HITS, Katz, e A column 
sums and PageRank with ct = 0.85. 





HITS 


Katz 


e cs 


PR 


6837 


6837 


6837 


6837 


2264 


6840 


6839 


6839 


6839 


8226 


6839 


6840 


6840 


6840 


8059 


6838 


6838 


6838 


6838 


8057 


6617 


6617 


6573 


6573 


4485 


6615 


6615 


6574 


6575 


5707 


6766 


6614 


6575 


6576 


8225 


6764 


6616 


6576 


6577 


6837 


6616 


6764 


6577 


6578 


6839 


6614 


6766 


6578 


6579 


6840 



observations is that the adjacency matrix A for this digraph is much closer to being 
symmetric than in the other cases. Indeed, the percentage of "bidirectional" edges in 
the wb-cs-stanford graph is 47.63%; the corresponding percentages for the abortion, 
computational complexity and death penalty graphs are just 2.72%, 2.97% and 4.02%, 
respectively. 

Interestingly, the (Reverse) PageRank results are now drastically different fron 
the ones provides by all the other measures in nearly all cases. The only (partial) 
exception is that PageRank finds nodes 6837, 6839 and 6840 to be among the top 10 
authorities; these three nodes are identified as the three most authoritative ones by 
the remaining methods. 

9. Approximating the matrix exponential. Several approaches are avail- 
able for computing the matrix exponential |26| . A commonly used scheme is the 
one based on Pade approximation combined with the scaling and squaring method 
[25| 126] . implemented in Matlab by the expm function. For an n X n matrix, this 
method requires 0(n 2 ) storage and 0(n 3 ) arithmetic operations; any sparsity in A, 
if present, is not exploited in currently available implementations. Evaluation of the 
matrix exponential based on diagonalization also requires 0(n 2 ) storage and 0(n 3 ) 
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operations. Furthermore, these methods cannot be easily adapted to the case where 
only selected entries (e.g., the diagonal ones) of the matrix exponential are of interest. 

For the purpose of ranking hubs and authorities in a directed network, only the 
main diagonal of e is required. This can be done without having to compute all the 
entries in e A . If some of the off-diagonal entries (communicabilities) are desired, for 
example those between the highest ranked hubs and/or authorities, it is also possible 
to compute them without having to compute the whole matrix e , which would be 
prohibitive even for a moderately large network. We further emphasize that in most 
applications one is not so much interested in computing an exact ranking of all the 
nodes in a digraph, but only in identifying the top k ranked nodes, where the integer 
k is small compared to n (for example, k — 10 or k = 20). It is highly desirable to 
develop methods that are capable of quickly identifying the top k hubs/ authorities 
without having to compute accurate hub/authority scores for each node. 

Efficient, accurate methods for estimating (or, in some cases, bounding) arbitrary 
entries in a matrix function f(A) have been developed by Golub, Mcurant and col- 
laborators (see |24j and references therein) and first applied to problems of network 
analysis by Benzi and Boito in [2]; see also [6]. Here we limit ourselves to a brief de- 
scription of these methods, referring the reader to [2] and [24] for further details. Let 
A be a real, symmetric, n x n matrix and let / be a function defined on the spectrum 
of A. Consider the eigendecompositions A = QAQ T and f(A) = Qf(A)Q T , where 
Q = . . . , cf) n ] and A = diag (Ai, . . . , A„); here we assume that the eigenvalues of A 
are ordered as Ai < . . . < A„. For given vectors u and v we have 

n 

u T f(A)v = u T Qf(A)Q T v = w T f(A)z = £ f(X k )w k z k , (9.1) 

fe=i 

where w = Q T u = (w k ) and z = Q T v = (zk)- In particular, for f(A) = e A we obtain 

n 

u T e A v = e Xk w k z k . (9.2) 
fe=i 

Choosing u = v = a (the vector with the ith entry equal to 1 and all the remaining 
ones equal to 0) we obtain an expression for the subgraph centrality of node i: 

n 

SC(i) : X '"•'•'"•<• 

fe=i 

where (j>k,i denotes the ith component of vector cf> k . Likewise, choosing u — e, and 
v = ej we obtain the following expression for the communicability between node i 
and node j: 

n 

C(i,j) :=^2e Xk <t)ka<pkj. 
fe=l 

Analogous expressions hold for other matrix functions, such as the resolvent. 

Hence, the problem is reduced to evaluating bilinear expressions of the form 
u T f(A)v. Such bilinear forms can be thought of as Ricmann- Sticltjcs integrals with 
respect to a (signed) spectral measure: 

,b ( 0, if A < a = Ai, 

u T f(A)v = / f(X)dfj,(X), fi(X) = < Y?k=i w k z k, if A 4 < A < X t+1 , 
Ja { ELi w kz k , if b = X n < A. 
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This integral can be approximated by means of a Gauss-type quadrature rule: 

f f(x)dfx(x) = cjfitj) + v *f^) + R ifi ( 9 - 3 ) 
Ja j=i fe=i 

where R[f] denotes the error. Here the nodes {tj}^ =1 and the weights {cj}*j =1 are 
unknown, whereas the nodes {rk}1 =1 are prescribed. We have 

• q = for the Gauss rule, 

• q = 1, t\ = a or t\ = b for the Gauss-Radau rule, 

• q = 2, T\ = a and T2 = b for the Gauss-Lobatto rule. 

For certain matrix functions, including the exponential and the resolvent, these 
quadrature rules can be used to obtain lower and upper bounds on the quantities of 
interest; prescribing additional quadrature nodes leads to tighter and tighter bounds, 
which (in exact arithmetic) converge monotonically to the true values [53] . The eval- 
uation of these quadrature rules is mathematically equivalent to the computation 
of orthogonal polynomials via a three-term recurrence, or, cquivalcntly, to the com- 
putation of entries and spectral information of a certain tridiagonal matrix via the 
Lanczos algorithm. Here we briefly recall how this can be done for the case of the 
Gauss quadrature rule, when we wish to estimate the ith diagonal entry of f{A). It 
follows from (|9.3p that the quantity of interest has the form X^=i c jf(tj)- This can 
be computed from the relation (Theorem 3.4 in 

p 

^Z^jf^j) = ef/(J P )ei, 

where 

/ wi 7i \ 

7l w 2 72 

7p-2 w p _x 7 P _i 
\ 7 P -i w p / 

is a tridiagonal matrix whose eigenvalues are the Gauss nodes, whereas the Gauss 
weights are given by the squares of the first entries of the normalized eigenvectors of 
J p . The entries of J p are computed using the Lanczos algorithm with starting vectors 
:r_i =0 and xq = e^. Note that it is not required to compute all the components of 
the eigenvectors of J p if one uses the Golub-Welsch QR algorithm; see [UJ . 

For small p (i.e., for a small number of Lanczos steps), computing the (1,1) entry 
of f(J p ) is inexpensive. The main cost in estimating one entry of f(A) with this 
approach is associated with the sparse matrix- vector multiplies in the Lanczos algo- 
rithm applied to the adjacency matrix A. If only a small, fixed number of iterations 
are performed for each diagonal element of f(A), as is usually the case, the compu- 
tational cost (per node) is at most 0(n) for a sparse graph, resulting in a total cost 
of 0(n 2 ) for computing the subgraph ccntrality of every node in the network. If only 
k < n subgraph centralities are wanted, with k independent of n, then the overall cost 
of the computation will be 0(n) provided that sparsity is carefully exploited in the 
Lanczos algorithm and that only a small number p of iterations (independent of n) is 
carried out. Note, however, that depending on the connectivity characteristics of the 
network under consideration, the prefactor in the 0[n) estimate could be large. The 
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Table 9.1 

The number of iterations necessary for the top 10 hubs or authorities to be determined (not 
necessarily in the correct order). 



Datasct 


hub (lower bound) 


hub (upper bound) 


Abortion 


> 40 


> 40 


Comp. Complex. 


3 


3 


Death Penalty 


5 


3 


Stanford 


8 


8 



Dataset 



authority (lower bound) 



authority (upper bound) 



Abortion 
Comp. Complex. 
Death Penalty 
Stanford 



2 
4 
4 
7 



algorithm can be implemented so that the storage requirements are 0{n) for a sparse 
network — that is, a network in which the total number of links grows linearly in the 
number n of nodes. 

When applying the approach based on Gauss quadrature rules to the 2n x 2n 
matrix A, only matrix-vector products with A and its transpose are required, just 
like in the HITS algorithm. If only the hub scores are wanted, it is also possible 
to apply the described techniques to the symmetric matrix AA T using the function 
/(A) = cosh(v'A); the same applies if only the authority scores are wanted, working 
this time with the matrix A T A. The problem with this approach is that only estimates 
(rather than increasingly accurate lower and upper bound) can be obtained, due to 
the fact that the function /(A) = cosh(-\/A) is not strictly completely monotonic on 
the positive real axis. We refer to [4] for details. In our experiments we always work 
with the matrix A, since we are interested in computing both hub and authority 
scores. 

9.1. Test results. Accurate evaluation of all the diagonal entries of e A using 
quadrature rules may be too expensive for truly large-scale graphs. In most applica- 
tions, fortunately, it is not necessary to rank all the nodes in the network: only the 
top few hubs and authorities are likely to be of interest. When using quadrature rules, 
the number of quadrature nodes (Lanczos iterations) required to correctly rank the 
nodes as hubs or authorities varies and depends on both the eigenvalues of e -4 and 
how close the diagonal entries are in value. If the rankings of the nodes are very close, 
it can take many iterations for the ordering to be exactly determined. However, since 
estimates for diagonal entries are calculated individually, once the top 10 (say) nodes 
have been identified, additional iterations can be performed only on these nodes in 
order to determine their exact ranking. 

Our approach exploits the monotonicity of the Gauss-Radau bounds: as soon as 
the lower bound for node i is above the upper bounds for other nodes, we know that 
node i will be ranked higher than those othe nodes. This observation leads to a simple 
algorithm for identifying the top-/c nodes. The number of Lanczos iterations per node 
necessary to identify the top k = 10 hubs and authorities, using Gauss-Radau lower 
and upper bounds, for the four data sets from section [8] is given in Table 19.11 Our 
implementation is based on Meurant's Matlab code [Mj, From the table it can be 
seen that, in most cases, only 2-5 iterations per node are needed. An exception is the 
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Table 9.2 

The number of iterations necessary for the top 10 hubs or authorities to be ranked in the top 30. 



Dataset 



Abortion 
Comp. Complex. 
Death Penalty 
Stanford 



hub (lower bound) 



2 
2 
7 



hub (upper bound) 



Dataset 



Abortion 
Comp. Complex. 
Death Penalty 
Stanford 



authority (lower bound) 



authority (upper bound) 



determination of the top 10 hubs of the abortion data set, for which the number of 
iterations is large (> 40). This is due to a cluster of nodes (nodes 960 and 968-990) 
that have nearly identical hub rankings. These nodes' scores agree to 15 significant 
digits. However, for most applications, if a subset of nodes are so closely ranked, their 
exact ordering may not be so important. Table 19.21 reports the number of Lanczos 
iterations neeeded for the top k = 10 hubs and authorities to be ranked at least 
in the top 30. Here, the number of iterations per node needed is never more than 
7. The total cost is thus 0{n) Lanczos iterations, again leading to an 0(n 2 ) overall 
complexity. Various enhancements can be used to reduce costs, including the use of 
sparse-sparse mat-vecs in the Lanczos iteration, and the exclusion of nodes with zero 
out-degree (for hub computations) and zero in-degree (for authority computations) 
from the top-k calculations. It is also safe to assume that in most cases of interest, 
one can also exclude nodes with in- and out-degree 1 from the computations, leading 
to further savings. 

10. Conclusions and outlook. In this paper we have presented a new approach 
to ranking hubs and authorities in directed networks using functions of matrices. Bi- 
partization is used to transform the original directed network into an undirected 
one with twice the number of nodes. The adjacency matrix of the bipartite graph 
is symmetric, and this allows the use of subgraph centrality (and communicability) 
measures for undirected networks. We showed that the diagonal entries of the matrix 
exponential provide hub and authority rankings, and we gave an interpretation for 
the off-diagonal entries (communicabilitics) . Unlike HITS, the results are indepen- 
dent of any starting vectors; and unlike the Katz-bascd ranking schemes, there is no 
dependency on an arbitrary parameter. 

Several examples, both synthetic and corresponding to real data sets, have been 
used to demonstrate the effectiveness of the proposed ranking algorithms relative 
to HITS and to other ranking schemes based on the matrix resolvent and on the 
exponential of the adjacency matrix of the original digraph. Our experiments indicate 
that our method results in rankings that are frequently different from those computed 
by HITS, at least in the absence of large gaps between the dominant singular value 
of the adjacency matrix A and the remaining ones. This is to be expected, since our 
method uses information from all the singular spectrum of the network, not just the 
dominant left and right singular pairs. 

As usual in this held, there is no simple way to compare different ranking schemes, 
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and therefore it is impossible to state with certainty that a ranking scheme will give 
"better" results than a different scheme in practice. It is, however, certainly the case 
that the method based on the exponential of A takes into account more spectral 
information than HITS does; moreover, the rankings so obtained are unambiguous, 
in that they do not depend on an the choice of an initial guess or on a tuneable 
parameter. As we saw, the latter is a weak spot of Katz-like methods, and a similar 
case can be made for PageRank and Reverse PagcRank. 

Compared to HITS, the new technique has a higher computational cost. We 
showed how Gaussian quadrature rules can be used to quickly identify the top ranked 
hubs and authorities for networks involving thousands of nodes. We note that such 
schemes require a symmetric input matrix and are not readily applicable to nonsym- 
metric matrices, since in this case one can only hope for estimates instead of lower 
and upper bounds. 

Future work should include a more efficient implementation and tests on larger 
networks. It is likely that the proposed approach based on Gaussian quadrature 
will prove to be too expensive for truly large-scale networks with millions of nodes. 
We hope to explore techniques similar to those presented in [5] and [33] in order to 
extend our methodology to truly large-scale networks. Another relevant question is 
the study of the rate of convergence of the Lanczos algorithm for estimating bilinear 
forms associated with adjacency matrices of graphs of different types. 
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